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ABSTRACT 

We carry out a set of self-consistent TV-body calculations to investigate how important 
the velocity anisotropy in non-spherical dark-matter halos is for dynamical friction. For 
this purpose we allow satellite galaxies to orbit within flattened and live dark-matter 
haloes (DMHs) and compare the resulting orbit evolution with a semi-analytic code. 
This code solves the equation of motion of the same satellite orbits with mass loss 
and assumes the same DMH, but either employs Chandrasekhar's dynamical friction 
formula, which does not incorporate the velocity anisotropy, or Binney's description of 
dynamical friction in anisotropic systems. In the numerical and the two semi-analytic 
models the satellites are given different initial orbital inclinations and orbital eccen- 
tricities, whereas the parent galaxy is composed of a DMH with aspect ratio qh — 0.6. 

We find that Binney's approach successfully describes the overall satellite de- 
cay and orbital inclination decrease for the whole set of orbits, with an averaged 
discrepancy of less than 4 per cent in orbital radius during the first 3 orbits. If Chan- 
drasekhar's expression is used instead, the discrepancy increases to 20 per cent. Bin- 
ney's treatment therefore appears to provide a significantly improved treatment of 
dynamical friction in anisotropic systems. 

The velocity anisotropy of the DMH velocity distribution function leads to a 
significant decrease with time of the inclination of non-polar satellite orbits. But at 
the same time it reduces the difference in decay times between polar and coplanar 
orbits evident in a flattened DMH when the anisotropic DMH velocity distribution 
function is not taken into account explicitly. Our TV— body calculations furthermore 
indicate that polar orbits survive about 1.6 times longer than coplanar orbits and 
that the orbital eccentricity e remains close to its initial value if satellites decay slowly 
towards the galaxy centre. However, orbits of rapidly decaying satellites modelled with 
the semi-analytic code show a strong orbital circularisation (e < 0) not present in the 
N-body computations. 

Key words: stellar dynamics - methods: N-body simulations- methods: analytical 
- galaxies: kinematics and dynamics - galaxies: haloes - galaxies: dwarf 



1 INTRODUCTION 

According to current ideas galaxy formation began with 
small-amplitude Gaussian fluctuations at the early stages 
of the Universe. In hierarchical cosmological models, these 
fluctuations decrease with increasing scales, resulting first in 
the formation of low-mass objects that may merge, building 
up ever more massive structures. The shape and morphol- 
ogy of these objects are strongly dependent on the cosmo- 
logical models, as one can conclude from N-body computa- 
tions, although none of them predict spherical structures. 
The most successful hierarchical theory is the Cold Dark 



Matter model (CDM). In this framework, aspherical bound 
dark matter haloes (DMHs) form as a result of gravitational 
clustering. The inclusion of gas dynamics in the CDM sim- 
ulations (Udry & Martinet 1994, Dubinsky 1994) results 
to a Gaussian distribution of DMH density aspect ratios, 
qh = c/a > 0, where c and a are the minor and major axes 
of an oblate spheroid, of mean < qn >= 1/2 and dispersion 
equal to 0.15 (Dubinsky 1994). The degree of asphericity 
depends on the dark matter nature. For example, numerical 
computations based on the ACDM predict halo axis-ratios 
of q h ~ 0.7 (Bullock 2001), Hot Dark Matter models lead to 
haloes as round as qh = 0.8 (Peebles 1993), whereas candi- 
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dates such as cold molecular gas (Pfenniger, Combes & Mar- 
tinet 1994) and massive decaying neutrinos (Sciama 1990) 
may produce halo profiles as flattened as qn = 0.2. 

Observationally, measuring galaxy axis-ratio is compli- 
cated: (i) Stellar kinematics: Oiling & Merrifield (2000) ob- 
tain an axis-ratio of qh « 0.8 for our Galaxy. This method 
has the disadvantage of having access to information of our 
Galaxy only on small scales, (ii) The flying gas layer method 
(Oiling 1996, Becquaert, Combes & Viallefond 1997): as- 
suming that the HI emission comes from gas in hydro- 
static equilibrium in galactic potentials axis-ratios as low as 
qh ~ 0.3 are obtained for the galaxies NGC 891 and 4244, 
(iii) Warping gas layer: Hofner & Sparke (1994) find axis- 
ratios of approximately 0.7 for NGC 2903 and of qn ~ 0.9 
for NGC 2841, 3198, 4565 and 4013, (iv) X-ray isophotes: 
Boute & Canizares (1998) measure values of qh « 0.5 for 
NGC 3923, 1332 and 720, (v) Polar ring galaxies: Arnaboldi 
ct al. (1993) and Sackett et al. (1994) find an axis-ratio 
of q h « 0.3 for NGC 4650A, 0.5 for the galaxy A0136- 
0801 and 0.6 for AM2020-504, (vi) Precessing dusty discs: 
Steinman-Cameron, Kormendy & Durisen (1992) measure 
an axis-ratio of 0.9 for the galaxy NGC 4753. 

Another method, which we focus on here, is the anal- 
ysis of satellite dynamics. There are two main different ap- 
proaches to infer the halo shape from satellites. First, one 
may attempt to reproduce the observed tidal streams of 
Milky Way satellites as done, for instance, by Ibata et al. 
(2001) and Law et al. (2003) who, from measurements of ve- 
locity, position and structure of the Sagittarius dwarf galaxy, 
constrain the initial parameter space and, subsequently, cal- 
culate in detail the satellite mass loss. They find that the 
Milky Way halo potential cannot be more flattened than 
« 0.8, otherwise tidal streams would be too spread out and 
thick compared to the observations due to orbital preces- 
sion. The second approach is a statistical study of satellite 
distribution around spiral galaxies. Holmberg (1969), Zarit- 
sky & Gonzalez (1999), Prada et al. (2003) and Sales & 
Lambas (2003) point out that satellites around disc galax- 
ies are found more often aligned with the poles of the host 
galaxy, the so-called 'Holmberg effect', whereas Quinn & 
Goodmann (1986) find in their iV-body study that the disc 
alone cannot account for the original statistical distribution 
of Holmberg's data. A remedy may be sought in the form 
of an extended non-spherical DMH. An anisotropic velocity 
(and mass) distribution will cause a satellite's orbit to align 
with the axes of the velocity ellipsoid of the host galaxy 
(Penarrubia, Kroupa & Boily 2001, hereinafter PKB). In 
both schemes, a large number of numerical calculations is 
needed. In the former approach one should integrate several 
initial orbital parameters to find the best fit to the observed 
satellite characteristics, whereas in the latter approach an 
initial satellite sample should statistically reproduce the dis- 
tributions expected from cosmological models of DMH for- 
mation. So far this is prohibitively time-expensive by means 
of any of the present N-body algorithms. 

Several studies of satellite decay have shown that, 
in spherical systems, Chandrasekhar's dynamical friction 
(Chandrasekhar 1943) is accurate enough if the Coulomb 
logarithm remains as a free parameter to fit to the N-body 
data (e.g. van den Bosch et al. 1999, Colpi et al. 1999) since 
it also depends on the code parameters and the number of 
particles employed. For instance, Prugniel & Combes (1992) 



and Whade & Donner (1996) find that dynamical friction is 
artificially increased due to numerical noise if the particle 
number is small. Semi-analytic methods that include Chan- 
drasekhar's dynamical friction have been demonstrated to 
reproduce accurately the overall evolution of satellite galax- 
ies (e.g Velazquez & White 1999, Taylor & Babul 2001) and, 
therefore, represent a useful tool in order to carry out ex- 
tensive studies with a large parameter space. 
It is, however, still unclear how the inhomogeneity of the 
system distribution affects the satellite dynamics. Whereas 
Del Popolo (2003), Del Popolo & Gambera (1998) and Maoz 
(1993) show that dynamical friction increases the steeper the 
density profile is, Just & Penarrubia (2002), following the 
theoretical scheme proposed by Binney (1977), hereinafter 
B77, find a negligible effect on the satellite orbit due to the 
symmetry of the inhomogeneous terms of dynamical friction. 

Although Chandrasekhar's formula, which assumes an 
isotropic velocity distribution, reproduces accurately dy- 
namical friction in spherical systems, an analytical study 
of Statler (1991) shows that in the case of Stackel potentials 
the velocity anisotropy produce strong effects on the satel- 
lite orbit. The N-body computations of PKB confirm that 
Chandrasekhar's formula cannot account for the resulting 
satellite decay and evolution if the halo is flattened. The aim 
of this paper is to implement a semi-analytic scheme capable 
of tracking the dynamical evolution of substructures within 
flattened as well as spherical DMH's. With this purpose in 
mind, we implement in our code the analytic expressions of 
dynamical friction in systems with anisotropic velocity dis- 
persions suggested by Binney (B77) and also used by Statler, 
which reproduces Chandrasekhar's for null anisotropy in ve- 
locity space. We as well compare the results of using Chan- 
drasekhar's formula in axi-symmetric systems to determine 
the effects of the velocity anisotropy on satellite decay. 

Section 2 introduces the models. In Section 3 we provide 
the code and galaxy parameters. We outline the dynamical 
friction approaches in Section 4 whereas in Section 6 we pro- 
pose a simple technique to fit a free parameter (in our case 
the Coulomb logarithm) to the N-body data. In Section 7 we 
study how flattened DMHs affect satellite decay and calcu- 
late the degree of accuracy of Binney's formula. The paper 
concludes with Section 8. 



2 GALAXY AND SATELLITE PARAMETERS 

Our DMH model is that used by PKB to facilitate an inter- 
comparison of the disc and bulge effects on satellite decay 
(Penarrubia 2003). Here we do not add the bulge and disc 
components in order to distill the dependency of dynamical 
friction on the velocity anisotropy in a spheroidal DMH. 

In order to minimise computational time when con- 
structing flattened DMHs, we apply a highly- efficient tech- 
nique using multi-pole potential expansions to tailor the 
local velocity ellipsoid to the required morphology (Boily, 
Kroupa & Penarrubia 2001). The MaGalie code scales lin- 
early with particle number and hence we can construct flat- 
tened DMHs consisting of ^ 10 6 particles or more, in a short 
computational time. 

Following PKB, the flattened DMH is described by a 
non-singular isothermal profile which, in cylindrical coordi- 
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nates (for m 2 [u = 0] see eq. 7), can be described as 
M h a exp[-m 2 (0)/r 2 ut ] 



Ph[m\Q)] 



2n*/ 2 r cut m 2 (0)+ 7 2 

with m 2 (0) = R 2 + z 2 /q 2 



(1) 



M h being the DMH mass, r cut the cut-off radius and 7 the 
core radius, and 



{l-0F/3exp(/3 2 )[l- 
1 + + (tt - 



erf^)]}"^ 
2)/3 2 + 0(/3 3 ) 



(2) 



where /3 = 7/r cut ^ 1/24 in our calculations. For /3 = 1/24 
we find a ~ 1.076 — + 1 already and hence thereafter we set 
q = 1 in our analysis. 

We use self-consistent King models (King 1966) to rep- 
resent our dwarf galaxies. These models fit early- type dwarf 
galaxies (Binggeli et al. 1984). For a comparison with the 
work of PKB we adopt c = log 10 (r t /r c ) = 0.8, where r c and 
r t are the core and tidal radii, respectively. 

To construct the models we choose the satellite mass 
M s and r t . The tidal radius is determined by computing 
the density contrast, p s (rt) /~Pg(r a ) ~ 3, at the apo-centric 
distance (r a = 55 kpc) at t = 0, Ps( r ) being the averaged 
density of the galaxy (same procedure as Velazquez & White 
1999). This guarantees that all satellite particles are bound 
at t = 0. 

We employ the system of units of PBK, which refers 
to the parameters of the Milky Way disc. Adopting Md = 
Rd = G = 1 and according to Bahcall, Smith & Soneira 
(1982), M d = 5.6 x 10 10 M o and R d = 3.5 kpc, the time and 
velocity units are, respectively, 1.3 x 10 7 yr and 262 kms -1 . 

The values of the galaxy and the satellite parameters 
can be found in Table 1. The parent galaxy corresponds to 
the model G2 of PKB, with the difference that we remove 
the disc and bulge components here. 



3 NUMERICAL CALCULATIONS 
3.1 Code parameters 

The numerical experiments were carried out by using the 
code Superbox which is a highly efficient particle mesh- 
algorithm based on a leap-frog scheme (for a detailed de- 
scription see Fellhauer et al. 2000). Superbox has already 
been implemented in extensive studies of satellite disrup- 
tion by Kroupa (1997), Klessen & Kroupa (1998) and PKB. 
The program calculates the accelerations using a high order 
NGP ('nearest grid point') force calculation scheme based 
on the second derivatives of the potential. A self-consistent 
system of several galaxies can be treated by forming sub- 
grids (3D 'boxes') which follow the motion of each galaxy. 
Each sub-grid has three levels of resolution, the two finest 
levels co-move with the galaxies allowing a high-resolution 
calculation of the forces acting on the particles, whereas the 
third one covers the local 'Universe'. The finest levels are 
centred on the density maximum of the galaxy, which is re- 
computed at every time-step. 

The code parameters are those of PKB. In that paper a 
detailed description of the system and the grid structure is 
presented, whereas here we merely give a brief description 
of the N-body parameters. Our integration time step is 0.39 





Symbol 


Valuc(ph.u) 


Value (m.u) 


DMH 


N h 


1 400 000 




(H2) 










M h 


7.84 x lO n M 


14.00 




Qh 


0.60 


0.60 




7 


3.5 kpc 


1.00 




»*cut 


84.00 kpc 


24.00 


Satellite 


N s 


40 000 




(SI) 










M s 


5.60 x 1O 9 M 


0.10 




*(o)K 2 


5.00 


5.00 




r c 


0.67 kpc 


0.19 




rt 


7.24 kpc 


2.07 




c 


1.03 


1.03 




< r > 


1.64 kpc 


0.47 




(JO 


eO.SOkms- 1 


0.23 



Table 1. Primary galaxy and satellite models. The DMH has an 
aspect ratio = 0.6. For the satellite model: 'f(O) = $(rt) — 
<1?(0), ^(O) are the central potential and <J>(ri) the potential at 
the tidal radius (following the notation of Binney & Tremaine 
1987); o"o is the velocity dispersion at the centre, and < r > 
the average radius of the satellite. The units are such that Ph.u. 
means 'physical units' and m.u. 'model units'. Nh and N s are the 
number of particles used to represent the DMH and the satellite, 
respectively. 



Myr which is about 1 /40th the dynamical time of our satel- 
lite. We have three resolution zones, each with 64 3 grid-cells: 
(i) The inner grid covers out to 3 radial halo scale- lengths, 
providing a resolution of 350 pc per grid-cell, (ii) The middle 
grid covers the whole galaxy, with an extension of 24 halo 
scale-lengths (84 kpc), giving a resolution of 2.8 kpc per 
grid-cell. The satellite always orbits within this grid except 
at the very late stages of its evolution, avoiding cross-border 
effects (see Fellhauer et al. 2000) and (iii) The outermost 
grid extends to 348 kpc and contains the local universe, at a 
resolution of 11.6 kpc. As for the satellite grid-structure, the 
resolutions are 816 pc per grid-cell for the inner grid that 
extends to 24.48 kpc, 1.2 kpc per grid-cell for the middle 
grid which extends to 36 kpc, and 11.6 kpc per grid-cell for 
the outermost grid that covers the local universe. 

The selection of grid parameters ensures the conserva- 
tion of energy and angular momentum for satellites in iso- 
lation over times as long as our calculations to better than 
1% for all the models. 



3.2 Orbital parameters 

The parameter space of satellite galaxies is extremely large. 
A complete survey should account for different satellite 
masses, apo-galacticon distances, orbital eccentricities and 
inclinations. ..etc. In this paper, we carry out a set of cal- 
culations selecting those parameters that best reflect the 
effects of the velocity dispersion anisotropy on satellite de- 
cay. These parameters are: (i) the initial orbital inclination 
(i), defined as the angle between the initial angular momen- 
tum vector of a satellite and the axis perpendicular to the 
axi-symmetry plane (selected as the z-axis). We expect or- 
bital inclination to decrease in time as predicted by Binney. 
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Name 


Gal. 


Sat. 


i 


e 


r p 


r a 




model 


model 






[kpc] 


[kpc] 


H2S100 


H2 


SI 


0° 


0.5 


18 


55 


H2S130 


H2 


SI 


30° 


0.5 


18 


55 


H2S145 


H2 


SI 


45° 


0.5 


18 


55 


H2S160 


H2 


SI 


60° 


0.5 


18 


55 


H2S190 


H2 


SI 


90° 


0.5 


18 


55 


H2S100c 


H2 


SI 


0° 


0.3 


30 


55 


H2S130c 


H2 


SI 


30° 


0.3 


30 


55 


H2S145c 


H2 


SI 


45° 


0.3 


30 


55 


H2S160c 


H2 


SI 


60° 


0.3 


30 


55 


H2S190c 


H2 


SI 


90° 


0.3 


30 


55 


H2S100e 


H2 


SI 


0° 


0.7 


10 


55 


H2S130c 


H2 


SI 


30° 


0.7 


10 


55 


H2S145e 


H2 


SI 


45° 


0.7 


10 


55 


H2S160e 


H2 


SI 


60° 


0.7 


10 


55 


H2S190c 


H2 


SI 


90° 


0.7 


10 


55 



Table 2. The numerical experiments. The peri- and apo-galactica 
are r p and r a , respectively, and e = (r a — r p )/(r a + r p ) is the 
orbital eccentricity. 



where (vr,v z ) are the coordinates of the satellite velocity in 
this frame. 

As Binney shows, a body with mass M s will suffer a 
decrease of its orbital inclination whenever B z > Br (oblate 
halo). If the orbit is either coplanar or polar, the inclination 
remains constant since, respectively, the perpendicular and 
the planar component of v is zero. It is straight-forward to 
show that this expression reproduces Chandrasekhar's for 
e v = 0, 

F ch = -4irGM s p h [m 2 (0)}lnA[erf(X) - ^e"* 2 ] ^, (4) 

where X = |v s |/-s/2cr. 

One important aspect to note is that both expressions 
of dynamical friction include here an anisotropic halo den- 
sity, denoted by ph[m 2 (0)] — ph[R 2 + z 2 /q\\ in cylindrical 
coordinates. This is the "local approximation" made here 
but we note that the derivation of eq. 4 assumes a uniform, 
infinitely extended background medium. In practice, the lo- 
cal approximation implies the only difference between both 
expressions to be the anisotropic terms of the velocity dis- 
tribution. 



We note that all the calculations proceed with the same 
orbital sense, but this is not important since the halo is 
non-rotating, (ii) The satellite's initial orbital eccentricity, 
denned as e = (r a — r p ) / (r a + r p ) , where r a , r p are the apo- 
and peri-galacticon, respectively. 

The parameters of the numerical experiments are listed 
in Table 2. 



4 HALO DYNAMICAL FRICTION 

Chandrasekhar's expression cannot explain some effects ob- 
served in N-body calculations of satellite decay within flat- 
tened haloes (PKB). Our aim is to check Binney 's approx- 
imation (B77) for systems with anisotropic velocity disper- 
sion. 

For simplicity, we reproduce here the analytic for- 
mulae employed throughout this study (for a detailed anal- 
ysis of the friction force see Penarrubia 2003 and Just & 
Penarrubia 2002). 

If the distribution function in velocity space is axi- 
symmetric, the zeroth order specific friction force is (B77) 



2V2^Ph[m 2 (0)]G 2 M s VT^ iJlnA D 
ti = 5 r>RVi, 



(3) 



2V2^Ph[m 2 {0)]G 2 M s VT=7 2 lnA 

r z = t, D X V Z 



where i — x,y and (<j r ,o z ) is the velocity dispersion el- 
lipsoid of a Schwarzschild distribution in cylindrical coordi- 
nates with constant ellipticity e 2 — 1 — (a z /aR) 2 , InA is the 
Coulomb logarithm of the halo and 



dq- 



exp( 



l + q 1 



(l + g) 2 (l-e? + g) 1 /2 ' 



B, 



f 

Jo 



dq- 



exp(- 



•'/2a 



i+q 



<' 2 < ) 



(l+ 9 )(l- e 2 +9 )3/2 



5 THE SEMI-ANALYTICAL CODE 

In order to analyse the accuracy of the analytic expressions 
of dynamical friction we have constructed a semi-analytic 
algorithm that solves the equation of motion of a point-mass 
satellite 



dt 2 



= F g + F df , 



(5) 



where F s is the specific force from the parent galaxy and 
Fdf that due to dynamical friction (eqs. 3 or 4). 

If the parent galaxy follows the fixed density profile of 
oq. (1), the specific force can be written in Cartesian coor- 
dinates (Chandrasekhar 1960) as 



F 9il = -27rGa; 



du 



(i +u ni + euw Ph[rn{u)] ' (6) 

roc 

F 9 , z = -2nGzj o {l+u){1 + el+u)3/2 PH\m 2 {u% 

where Xi = x,y, the ellipticity of the galaxy is e 2 h = 1 — q\, 
and 



2/ \ R 
m (u) = 



1+u 



1-ei + u 



(7) 



The algorithm employed to solve eq. (5) is based on 
the Bulirsch-Stoer method (for a complete description see 
Press et al. 1986), which provides high-accurate solutions 
with minimal computational effort. This method is based on 
an adaptive step-size scheme, thus, being ideal for systems 
with non-smooth potentials, as may be the case for satellites 
on highly eccentric orbits when disc and bulge components 
are included. 

For calculating Fdf (eqs. 3 and 4) the Coulomb log- 
arithm is treated as a free parameter to be fitted to the 
numerical orbits. The satellite mass M s (t) is obtained from 
the N-body data (see Section 7.1) and is treated as an input 
function. Numerical tests of orbits in a Keplerian potential 
show that energy and angular momentum are conserved to 
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10 -8 and 10 -9 per orbit, respectively, after choosing the ac- 
curacy of the Bulirsch-Stoer algorithm to be EPS= 1CP 5 . 
This value remains fixed for the calculations presented in 
this paper. An extensive description of our semi-analytic 
code can be found in Penarrubia (2003). 



6 DETERMINING THE COULOMB 
LOGARITHM 

The Coulomb logarithm is usually fit to numerical data with 
the aim of reproducing the overall orbital evolution by means 
of semi-analytic algorithms (e.g. Fellhauer et al. 2000). This 
procedure may actually be considered as a "calibration" of 
the semi-analytic code, which must be done carefully if a 
detailed inter-comparison between different schemes of dy- 
namical friction is desired. For that reason we present in 
what follows a method to describe the accuracy of the semi- 
analytic scheme. 

One possible way to quantify similarity of two orbits is 
via the quantity 



i=l 



i) 2 + a (ro)(U — U, n ) 



(8) 



r being the satellite position vector at the peri- and apo- 
galactica and t the time at which the satellite passes by 
these points. The subindex n denotes the numerical values 
and a(ro) the velocity dispersion of the DMH at the initial 
galacto-centre distance. The sum is over a given number of 
orbits k. 

The definition of \ measures the divergence of the nu- 
merical and semi-analytical satellite position. The term a At 
allows for the possibility that the orbital periods differ dur- 
ing the evolution, which would lead to a secular deviation . 
By definition, \ is equivalent to the discrepancy between the 
numerical and semi-analytical position evolution per orbit. 
The selection of the maximum and minimum galacto-centre 
distances for comparison permits a direct control over the 
orbital eccentricity evolution, although the measure of \ can 
be extended to other points without loss of generality. 

The value of k depends on the objectives of the study. 
For instance, if the aim is to find the best calibration for a 
large period of time, as it may be to reproduce the satellite 
decay in spiral galaxies, the fc-value must include a number 
of periods large enough, so that the overall evolution of sev- 
eral satellite orbits can be reproduced with a single value 
of In A. In this paper, we limit our fit to the first satellite 
periods, namely, k — 3, 4, for which the differences between 
both approaches of dynamical friction can be clearly seen. 

In Fig. 1 we plot \ f° r some of the experiments, con- 
cretely, those with inclinations 60°, 45° and 30° (rows), with 
eccentricities 0.3, 0.5 and 0.7 (columns). For each model, the 
semi-analytic code is employed to generate the satellite or- 
bit using Chandrasekhar's (dashed and dotted-dashed lines) 
and Binney's (full and dotted lines) formula to reproduce 
dynamical friction. This figure shows that Chandrasekhar's 
formula poorly describes the dependence of the satellite or- 
bit with the initial inclination, leading to a wider dispersion 
of the Coulomb logarithm values (for this range of inclina- 
tions, between 30° and 60°, InA G [0.9,2.8]). If Binney's 
expression is used, the variation of InA is highly reduced 
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Figure 1. The parameter \ for different orbital eccentricities 
and inclinations. Dotted lines denote the first 4 orbits, whereas 
solid lines the first 3 orbits, in both cases using Binney's expres- 
sions for dynamical friction. Dashed and dot-dashed lines repre- 
sent the results using Chandrasekhar's expression for k = 4 and 
3 orbits, respectively. The \ values are normalised to the initial 
apo-galacticon distance ro = 55 kpc. 



(InA £ [2.3,2.5]), which shows that this scheme provides a 
much better description of the effects of anisotropic velocity 
dispersion on satellite decay, independently of the orbital 
inclination. The variation range of the Coulomb logarithm 
becomes larger for a wider inclination spread, but x barely 
presents a dependence on the satellite eccentricity if Bin- 
ney's formula is applied. 

If Chandrasekhar's approach is used, the Coulomb log- 
arithm that leads to the best fit becomes smaller as the in- 
clination increases. Since dynamical friction is proportional 
to In A, the use of the average value implies an overestimate 
of the force for low inclinations and vice versa. The final 
average, 



(x/ro) = 



AW 



^Xi/ro, 



(9) 



over the iVi n A numerical experiments of Table 2 is plotted in 
Fig. 2. This figure shows the large discrepancies produced 
by Chandrasekhar's expression if the fit is for a large range 
of orbital inclinations and eccentricities, as expected. The 
minima of the curves determine the values of InA that lead 
to the best fits, which we summarise in Table 3. The values 
of Xmin denote the average error during the first k orbits 
associated with the fit. 

For the following analysis the value of the Coulomb 
logarithm implemented in our semi-analytic code will be 
found in Table 3. Looking at Fig. 1, we expect that Chan- 
drasekhar's friction will lead to more accurate fits for low 
(i ~ 30°) than for high inclined orbits after fixing InA = 2.2. 
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Figure 2. Average of the fitting parameters over the calculations 
of Table 2. 



Friction 


k 


InA 


Xmin/m 


Xmin/0 r ) 


Binney 


3 


2.4 


0.024 


0.0080 




4 


2.4 


0.038 


0.0095 


Chandrasekhar 


3 


2.1 


0.147 


0.049 




4 


2.2 


0.193 


0.048 



Table 3. Results of the fitting procedure applied to the numerical 
calculations of Table 2 for both formulae of dynamical friction. 
The fifth column gives relative deviation per orbit. 



Gm 3 



+ e 2 



(10) 



the softening being e ~ 0.23 mu = 0.8 kpc, which is the 
resolution of the inner grid focused on the satellite centre- 
of-density r s , and $ 9 the galaxy potential at this point, ne- 
glecting the tidal terms. Thus, all the particles of a satellite 
are assumed to feel the same external potential, which is a 
useful and sufficiently accurate approximation, taking into 
account that most of the bound particles are located very 
close to this point. 

Particles with Ei > are removed and the procedure 
repeated until only negative energy particles are left. As 
Johnston et al (1999) show, the energy criterium permits one 
to distinguish those particles that, though unbound, remain 
in orbits inside the tidal radius, which will escape from the 
satellite after some orbital periods. 

The mass is calculated each At = 0.312 Gyr, so that 
the semi-analytic code interpolates the value for intermedi- 
ate points at each time-step. The error is of the order of 
AM(t)/ At, going linearly with the mass loss. This means 
that the interpolation might introduce not negligible differ- 
ences at times where the mass loss is significant (i.e late 
times of the satellite evolution). In this study we are not 
concerned in detail with the late phases of evolution. 

In the right columns of Fig. 3a, b and c we plot the 
mass evolution for different orbits. Most of the satellites 
reach the inner most regions of the parent galaxy with a 
substantial fraction of their initial mass. A comparison of 
these curves with those of PKB (where a bulge and a disc 
component were included) shows that the baryonic subsys- 
tems of a galaxy induce a larger mass loss through tidal 
heating (e.g, Taylor & Babul 2001, Penarrubia 2003). PKB 
observe in their numerical experiments that all satellites 
with M s = O.lMd and ro = 55 kpc are destroyed before 
the remaining bound part of the satellite reaches the central 
region of the galaxy. 



7 THE VELOCITY ANISOTROPY EFFECTS 

In this Section, we discuss in more detail various aspects of 
satellite evolution. 



7.1 Satellite mass loss 

Tidal forces induce satellite mass loss. The satellite mass 
plays an important role in determining the ultimate fate 
of its evolution and survival (e.g. Velazquez & White 1999, 
Klessen & Kroupa 1998, Johnston et al. 1999). The value 
of M s (t) is treated as an input by the semi-analytical code, 
and is calculated using the self-consistent Superbox code. 

The mass remaining bound to the satellite, M„(t), is 
determined in the Superbox calculations by computing the 
potential energy <3>i < of each satellite particle presumed 
bound to the satellite, and its kinetic energy (Ti) in the 
satellite frame. Following PKB, particles with Ei = T + 
rris{&i + $ext) > are labelled unbound, where m 3 is the 
mass of one satellite particle (all have the same mass) and 



7.2 Satellite decay 

One of the most important effects of dynamical friction is 
the monotonic reduction of the orbital angular momentum 
and energy during the satellite's evolution that leads to a 
progressive decrease of the averaged galacto-centre distance. 
The computations carried out by PKB show a strong depen- 
dence of the decay time on the initial inclination that must 
be compared to analytic estimates. 

In Fig. 3a, b and c we plot the radius evolution (left 
columns) for those models with initial e = 0.5, 0.7 and 0.3, 
respectively. From this figure, we conclude that Binney's ex- 
pression clearly produces more accurate results than Chan- 
drasekhar's for the whole range of orbital inclinations. This 
result is not surprising due to the small dependence of the 
Coulomb logarithm on the inclination and eccentricity as 
shown in Fig. 1. Additionally, the value of InA that pro- 
duces the best fit for the first few orbits also succeeds in 
reproducing the decay time of the satellite. 

PKB observe that coplanar satellites suffer larger fric- 
tion than those following polar orbits, leading to survival 
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Figure 3. a: Radius and mass evolution for the models of Table 2 
with initial e = 0.5. Open circles denote the numerical evolution, 
whereas full and dotted lines represent the data obtained from the 
semi-analytic code using, respectively, Binney's (In A = 2.4) and 
Chandrasekhar's (In A = 2.2) expressions to reproduce dynamical 
friction. 



times 70% shorter. Due to the presence of a disc in their 
galaxy model, the contribution of the disc anisotropy on the 
decay differentiation as a function of the inclination can- 
not be directly measured. The calculations presented here 
(where the disc and bulge are removed) show survival times 
that range from 3.7 Gyr (coplanar orbits) to 6 Gyr (polar 
orbits) , using the same orbital parameters and halo flatten- 
ing as PKB. This implies a decay time difference of around 
60% between polar and coplanar satellites, which indicates 
that the disc contribution might be of the order of 10%. 
The effects of the disc on the satellite orbit can be found 
in Penarrubia (2003) and will be addressed in a following 
paper. 

Depending on the symmetry of the halo distribution, 
one can observe the following effects: 

• Spherical mass distribution and isotropic veloc- 
ity distribution: Satellites orbiting systems with a spheri- 
cal distribution function move on orbits that do not depend 
on their orientation with respect to the symmetry axis. 

• Flattened mass distribution and isotropic veloc- 
ity distribution: By means of the local approximation the 
value of dynamical friction is determined by the properties 
of the halo at the satellite's position, being reproduced by 
Chandrasekhar's formula. Our results indicate that the spa- 
tial asphericity leads to a strong differentiation of the satel- 
lite decay as a function of the orbital inclination. 

• Flattened mass distribution and anisotropic ve- 



Figure 3 - continued b: As Fig. 3a for those satellites of Table 2 
with initial e = 0.7. (Note that the time-axis has changed scale.) 
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Figure 3 - continued c: As Fig. 3a for those satellites of Table 2 
with initial e = 0.3. (Note that the time-axis has changed scale.) 
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locity distribution: The main influence of the velocity 
anisotropy on the satellite orbit is the secular reduction of 
the orbital inclination (see Subsection 7.3) and the reduc- 
tion of the spatial anisotropy effects, which is equivalent 
to Br < B z in Binney's formula (eq. 3, oblate systems). 
Taking into account the velocity anisotropy explicitly thus 
reduces the difference in orbital decay times between satel- 
lites on polar and coplanar orbits. As Fig. 3 shows, assum- 
ing an isotropic distribution in velocity space (Br — B z ) or, 
equivalcntly, using Chandrasekhar's formula to reproduce 
dynamical friction, leads to an overestimation of dynamical 
friction for low inclined satellites and to an underestimation 
for those satellites following highly inclined orbits. This is 
due to the fact that the anisotropy leads to a reduced effec- 
tive X for high inclined orbits yielding an enhanced friction 
force and vice versa for low inclination. This competes with 
the density effect due to the flattening, because highly in- 
clined orbits have lower mean density along the orbit com- 
pared to lower inclination. 

The evolution of the orbital radius of satellites with 
initial e = 0.3, 0.7 is plotted in Fig. 3b and c. As we can 
observe, Binney's approximation reproduces accurately the 
overall radius evolution independent of the initial eccentric- 
ity and orbital inclination. 




2 4 
t(Gyr) 



7.3 Evolution of the orbital inclination and 
eccentricity 

Orbits in non-spherical systems have inclinations (i = 
arccos[L z /L]) that do not remain constant but suffer pe- 
riodical oscillations due to nutation. In addition, the orbital 
plane precesses at a constant rate. Once the initial condi- 
tions are fixed, the amplitude and frequency of nutation 
and the precession rate remain constant if the friction force 
is removed from the equations of motion, whereas, if im- 
plemented, nutation and precession vary according to the 
angular momentum and radial distance evolution. Our in- 
terest focuses now on the effects induced by the velocity 
anisotropy on the satellite inclination during the orbit. 

Binney (B77) predicts a progressive reduction of i due 
to dynamical friction if the velocity dispersion ellipsoid is 
axi-symmetric (<jr,o z ) and o~r > a z . By symmetry the 
inclination decrease will not occur if the orbits are either 
coplanar (i = 0°) or polar (i — 90°). 

The inclination evolution of models with e = 0.5 is plot- 
ted in Fig. 4 (left column), where dotted lines denote the 
numerical data and solid and dashed lines the semi-analytic 
evolution if dynamical friction is modelled by Binney's and 
Chandrasekhar's formulae, respectively. This Figure shows 
the reduction of the mean value of i predicted by Binney and 
observed by PKB in their numerical experiments. After the 
satellite has sunk to the inner most region of the halo, the in- 
clinations are as low as 10° — 20° independent of their initial 
value. This large decrease of i is well reproduced by Binney's 
expression, although the nutation process shows discrepan- 
cies with the numerical result, which is connected with the 
not exact reproduction of the orbit. Despite the accurate de- 
scription of the overall decay process, this orbital mismatch 
is also observed when applying Chandrasekhar's expression 
in spherical systems (see Just & Penarrubia 2002). The fig- 



Figure 4. a: Inclination and eccentricity evolution for the models 
of Table 2 with e = 0.5. Left column: Dotted lines represent 
the N-body inclination evolution, full and dashed lines denote 
the use of Binney's and Chandrasekhar's equations, respectively. 
Right column: Numerical (dotted lines) against semi-analytical 
eccentricity evolution. Solid circles and open squares denote the 
implementation of Binney's and Chandrasekhar's expressions in 
the semi-analytic code, respectively. We note that, for clarity, the 
semi-analytic values of e are only plotted at the apo-centres. 



ure also confirms that the orbital inclination of coplanar and 
polar satellites remains constant. 

If dynamical friction is modelled by Chandrasekhar's 
formula, i.e the velocity distribution is assumed to be 
isotropic, the averaged value of i does not change during the 
orbit, which contradicts the numerical results. Thus, while 
the difference between the survival times of polar and copla- 
nar orbits is larger for flattened DMHs treated with Chan- 
drasekhar's dynamical friction formula (Fig. 3) , taking into 
account the DMH anisotropic velocity distribution function 
explicitly via Binney's formulae leads to an increase with 
time of the anisotropy of the satellite distribution due to 
the kinematical coupling of the satellites to the DMH ve- 
locity field. This clearly agrees with the fully self-consistent 
N— body computations reported here. 

In Fig. 4b and 4c (left columns) we plot the compar- 
ison for models with e = 0.7, 0.3, respectively. The results 
show barely a dependence on the eccentricity. It is interest- 
ing to note that, independently of e, orbits that are neither 
coplanar nor polar present large drops of the mean value 
of i. After the satellite sinks to the centre, the final orbital 
inclination lies in between 10-20° for all the models. 

We must emphasise the accuracy of Binney's formula 
in describing correctly the inclination decrease that satel- 
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Figure 4 - continued b: As Fig. 4a for models with e = 0.7. Note 
that the time-scale has a different value. 
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Figure 4 - continued c: As Fig. 4a for models with initially 
e = 0.3. Note that the time-scale has been rescaled. 



lites suffer in oblate systems. This is crucial for simulating 
properly satellite motions and for investigating satellite dis- 
tributions around spiral galaxies. 

Like the orbital inclination, the eccentricity is one of the 
orbital parameters that can be indirectly measured from ob- 
servations to determine a satellite's motion around a galaxy. 
In the right columns of Fig. 4 we show the comparison of the 
numerical eccentricity evolution with both semi-analytic ap- 
proaches. The analytic formula; of dynamical friction lead to 
a larger eccentricity decrease, which occurs mostly at late- 
times of the orbit, the so-called orbital circularisation, and 
becomes stronger for low inclined orbits, those that suffer 
higher dynamical friction. Fig. 4b and Fig. 4c indicate that 
the circularisation is more pronounced if the initial orbital 
eccentricity is higher and decreases for more circular orbits. 
Both dynamical friction expressions lead to a similar eccen- 
tricity evolution for the first few orbital periods, however at 
late-times the eccentricity exhibits a reduction not present 
in the numerical calculations that can be as high as ~ 80% 
for low inclined satellites on highly eccentric orbits (Fig. 4b, 
model H2S100e). 

7.4 Energy and angular momentum evolution 

A flattened system possesses two analytic constants of mo- 
tion, the energy and the component of the angular momen- 
tum perpendicular to the axi-symmetry plane (which we de- 
note as L z ). The total angular momentum L 2 — L\ + L 2 . 
is, however, not constant during a satellite orbit (see e.g 
Binney & Tremaine 1987), but has periodic variations that 
correspond to a precession and nutation of the orbital plane 
around the z-axis. 

Since the dynamical friction force has an opposite sense 
with respect to the satellite velocity, it decreases the angu- 
lar momentum and the energy which induces a monotonic 
sinking into the inner regions of the halo potential. The re- 
duction of angular momentum, therefore, implies an increase 
of the binding energy (in absolute value), since the potential 
grows with decreasing radius. Due to the small magnitude 
of dynamical friction compared to the mean field force, we 
expect an easier comparison between numerical and semi- 
analytic data by the slow variation of L z and E during the 
orbit. 

In Fig. 5 we plot the changes of specific E and L z due 
to dynamical friction for the models with e = 0.5. The re- 
sults are equivalent to those of the radial evolution. Due to 
our selection of In A as the average over those Coulomb log- 
arithms that lead to the best fit for each particular model, 
Chandrasekhar's formula overestimates dynamical friction 
for low inclined orbits and underestimates it for highly in- 
clined orbits. For orbits with i < 30°, this appears as a 
stronger reduction of the z-component of angular momen- 
tum and, equivalently, a large increase of the binding energy. 
The effect is contrary for satellites with i > 30°. 

This figure illustrates how the kinetic energy of the 
satellite is lost via friction, being taken-up partially by halo 
particles and also being deposited in unbound satellite par- 
ticles. At the end of the simulation the angular momentum 
has a null value, i.e the satellite remains in the inner most 
part of the galaxy. 

It is interesting to note that the numerical evolution of 
energy and angular momentum presents small oscillations 
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Figure 5. Specific energy and angular momentum evolution dur- 
ing the orbits with e = 0.5. The numerical evolution is denoted 
by dotted lines, whereas the semi-analytic data is represented by 
solid and dashed lines if dynamical friction is modelled by Bin- 
ney's and Chandrasckhar's formulae, respectively. The quantities 
E and L z are normalised to the initial value. Note that for the 
case i = 90° one has L z = 0. 

during the orbit. This behaviour is due to the self-response 
of the halo to the satellite motion. Since SUPERBOX preserves 
the total energy and angular momentum, the halo also moves 
around the centre-of-mass of the system. Due to the com- 
plexity of the feedback, it cannot be reproduced analytically 
(the halo centre-of-mass is fixed at the coordinate origin in 
the semi-analytic code). 



8 CONCLUSIONS 

To asses the accuracy of Binney's equations (B77) and in 
order to reproduce the decay of satellites in flattened DMHs 
with a semi-analytical code, we have performed a set of 
self-consistent numerical experiments for different orbital 
inclinations and eccentricities of the satellite. The semi- 
analytic code incorporates dynamical friction either in terms 
of Chandrasekhar's expression or as Binney's formulae. Both 
treatments include the aspherical density profile by means of 
the local approximation. This means that the differences on 
the satellite motion induced by each treatment of dynamical 
friction comes from the anisotropy in velocity space, which 
is implemented in the analysis of B77. 

The accuracy of Binney's and Chandrasekhar's for- 
mulae in comparison with the numerical orbits is determined 
by the parameter X 2 = w XX^ 1-2 + a ' 2 ^-t 2 ) at the peri- and 



apo-centres for a given number k of orbits. If dynamical fric- 
tion is modelled by Binney's equation, this quantity shows 
discrepancies of approximately Xmin = 0.009ro per orbit af- 
ter averaging over the set of experiments and for the first 
three orbits, while Chandrasekhar's formula produces values 
of around x = 0.05ro per orbit (Table 3). 

We conclude that Binney's expression faithfully re- 
produces the process of dynamical friction in anisotropic 
systems. The fit is as accurate as that employing Chan- 
drasekhar's formula in isotropic systems (see Just & 
Penarrubia 2002). 

The comparison of orbits resulting from Chan- 
drasekhar's and Binney's expression of dynamical friction 
gives us the possibility to asses the effects of the DMH ve- 
locity anisotropy on satellite dynamics. We have demon- 
strated that, (i) if the density profile is in both equations 
p = p(R,z), where R, z are the cylindrical coordinates of 
the satellite position vector, the orbits generated by Chan- 
drasekhar's formula overestimate the decay time for polar 
orbits and underestimate it for coplanar ones for an over- 
all best-fit Coulomb logarithm (InA = 2.2). One effect of 
the velocity anisotropy is then to reduce the interval of de- 
cay times as a function of the orbital inclination. Binney's 
expression has been shown to reproduce accurately the nu- 
merical results independently of the initial eccentricity and 
with InA = 2.4. (ii) Dynamical friction in systems with 
anisotropic velocity distribution leads to a marked decrease 
of the orbital inclination i which is well reproduced by Bin- 
ney's expression. After the satellite sinks to the inner most 
region of the galaxy, i lies within 10-20°, independently of 
the initial value unless i ~ 90° (polar orbits), (iii) The en- 
ergy and angular momentum evolution as a function of the 
orbital inclination confirm the results of (i) and (ii). 

The semi-analytic eccentricity evolution, either employ- 
ing Chandrasekhar's or Binney's formula, shows the so- 
called circularisation process, defined as the progressive re- 
duction of e during the orbit. This variation is stronger 
for increasing friction (as for coplanar orbits or during the 
late-times of the evolution) and barely takes place in the 
self-consistent N— body calculations. The small A— body 
circularisation agrees with the results of van den Bosch 
et al. (1999). A possible solution may be sought in the 
position-dependence of the Coulomb logarithm, as proposed 
by Hashimoto, Funato & Makino (2003). However, despite 
the improvement in the description of the orbit at early- 
times, the scheme of Hashimoto et al. overestimates the 
satellite decay time for all the experiments (see Just & 
Penarrubia 2002). We must conclude that the reason for 
circularisation in the semi-analytic orbits is not yet fully 
understood. 
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